/********************************************************************
*  sample.cu
*  This is a example of the CUDA program.
*********************************************************************/
#include <treeCuda.h>

#define PI 3.14159265

__device__ SimpleVertex float_to_vertex(float3 pos, float3 nor, float2 uv)
{
	SimpleVertex out;
	out.Pos.x=pos.x;
	out.Pos.y=pos.y;
	out.Pos.z=pos.z;

	out.Nor.x=nor.x;
	out.Nor.y=nor.y;
	out.Nor.z=nor.z;

	out.Tex.x=uv.x;
	out.Tex.y=uv.y;

	return out;
}

__device__ float3 cross(float3 a, float3 b)
{
	float3 c=make_float3(a.y*b.z-a.z*b.y,
						 a.z*b.x-a.x*b.z,
 						 a.x*b.y-a.y*b.x);
	
	return c;
}

__device__ float dot(float3 a, float3 b)
{
	float c=a.y*b.y+a.x*b.x+a.z*b.z;
					
	
	return c;
}

__device__ void add(float3*t ,float3 a, float3 b)
{
	t->x=a.x+b.x;
	t->y=a.y+b.y;
 	t->z=a.z+b.z;

}

__device__ float3 sub(float3 a, float3 b)
{
	float3 c=make_float3(a.x-b.x,
						 a.y-b.y,
 						 a.z-b.z);
	return c;
}

__device__ void mul(float3* t,float3 a, float b)
{
	t->x=a.x*b;
	t->y=a.y*b;
 	t->z=a.z*b;
}

__device__ void rotate(float3 * t,float3 vec, float3 u, float angle)
{
		float a=u.x*u.x + (1 - u.x*u.x)*cos(angle);
		float b=u.x*u.y*(1-cos(angle)) - u.z * sin(angle);
		float c=u.x*u.z*(1-cos(angle)) + u.y * sin(angle);
		float d=u.x*u.y*(1-cos(angle)) + u.z * sin(angle);
		float e=u.y*u.y + (1 - u.y*u.y)*cos(angle);
		float f=u.y*u.z*(1-cos(angle)) - u.x * sin(angle);
		float g=u.x*u.z*(1-cos(angle)) - u.y * sin(angle);
		float h=u.y*u.z*(1-cos(angle)) + u.x * sin(angle);
		float i=u.z*u.z + (1 - u.z*u.z)*cos(angle);
	
	
	
		t->x=vec.x*a+ vec.y*d+ vec.z*g;
		t->y=vec.x*b+ vec.y*e+ vec.z*h;
 		t->z=vec.x*c+ vec.y*f+ vec.z*i;



}

__device__ void normalize (float3 * a)
{
	float l = sqrt(a->x*a->x + a->y*a->y + a->z*a->z);
	a->x/=l;
	a->y/=l;
	a->z/=l;
}


/************************************************************************/
/* Init CUDA                                                            */
/************************************************************************/
#if __DEVICE_EMULATION__

bool InitCUDA(void){return true;}

#else
bool InitCUDA()
{
	int count = 0;
	int i = 0;

	cudaGetDeviceCount(&count);


	for(i = 0; i < count; i++) {
		cudaDeviceProp prop;
		if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
			if(prop.major >= 1) {
				break;
			}
		}
	}

	cudaSetDevice(i);


	return true;
}

#endif





__global__ void generate(SimpleVertex *vertices, DWORD* indices, branch * b, int w, float r)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
  
	
	
	branch current_b=b[idx];
	int offset_v=idx * w * VPW; 
	int offset_i=idx * w * IPW;

	float angle =(2 * PI) / w;
	
	float3 direction= make_float3(current_b.direction[0], current_b.direction[1], current_b.direction[2]);
	float3 prevDirection=make_float3(current_b.prevDirection[0], current_b.prevDirection[1], current_b.prevDirection[2]);
	float3 position= make_float3(current_b.position[0], current_b.position[1], current_b.position[2]);
	normalize(&direction);


	float3 side;
	if(acosf(dot(direction, make_float3(1,0,0)))>0.1f)
		side=cross(direction, make_float3(1,0,0));
	else
		side=cross(direction, make_float3(0,1,0));
	normalize(&side);
	
	float3 endUp=cross(direction, side);
	normalize(&endUp);

	if(acosf(dot(prevDirection, make_float3(1,0,0)))>0.1f)
		side=cross(prevDirection, make_float3(1,0,0));
	else
		side=cross(prevDirection, make_float3(0,1,0));
	normalize(&side);

	float3 startUp=cross(prevDirection, side);
	normalize(&startUp);

	
	for(int  i=0; i<w; i++)
	{
		
		//Generate vertices and indices
		float3 u;
		rotate(&u, endUp, direction, i*angle);
		normalize(&u);

		float3 v;
		rotate(&v, startUp, prevDirection, i*angle);
		normalize(&v);
		

		float r_s=r;
		
		float3 start_radius;
		mul(&start_radius,v,(expf((float)current_b.level/5.0f)-1.0f)*r_s);
		
		float3 end_radius;
		mul(&end_radius,u,(expf((float)(current_b.level-1)/5.0f)-1.0f)*r_s);
		
		float3 length;
		mul(&length,direction, current_b.length);
		
		float3 a;
		add(&a,position,start_radius);
		float3 b;
		add(&b,position,end_radius);
		add(&b,b,length);

		SimpleVertex start;
		start.Pos=D3DXVECTOR3(a.x, a.y, a.z);
		start.Nor=D3DXVECTOR3(v.x, v.y, v.z);
		start.Tex=D3DXVECTOR2(i%2,0);
		start.Tan=D3DXVECTOR3(prevDirection.x,prevDirection.y,prevDirection.z);
		float3 t=cross(v,prevDirection);
		start.Bin=D3DXVECTOR3(t.x, t.y, t.z);

		SimpleVertex end;
		end.Pos=D3DXVECTOR3(b.x, b.y, b.z);
		end.Nor=D3DXVECTOR3(u.x, u.y, u.z);
		end.Tex=D3DXVECTOR2(i%2,1);
		end.Tan=D3DXVECTOR3(direction.x,direction.y,direction.z);
		t=cross(u,direction);
		end.Bin=D3DXVECTOR3(t.x, t.y, t.z);

		vertices[offset_v+(i*VPW)]=start;
		vertices[offset_v+(i*VPW)+1]=end;
		
		indices[offset_i+i*IPW]=  offset_v+ i*VPW;
		indices[offset_i+i*IPW+1]=offset_v+ i*VPW+1;
		indices[offset_i+i*IPW+2]=offset_v+ ((i+1)%w) *VPW + 1;
		indices[offset_i+i*IPW+3]=offset_v+ i*VPW;
		indices[offset_i+i*IPW+4]=offset_v+ ((i+1)%w) *VPW + 1;
		indices[offset_i+i*IPW+5]=offset_v+ ((i+1)%w) *VPW;
		
		
	}


	
}


__device__ void generateLeaf(SimpleVertex *vertices, int offset,  float3 pos, float3 direction, float3 up, float size)
{
	normalize(&up);
	normalize(&direction);

	float3 left;
	float3 right;

	rotate(&left, direction, up, (PI / 3));
	rotate(&right, direction, up, -(PI / 3));

	normalize(&left);
	normalize(&right);

	mul(&direction,direction, size);
	mul(&left, left, size);
	mul(&right, right, size);

	float3 c, v1, v2, v3, v4, v5;

	add(&c, pos, direction);

	add(&v1, pos, left);

	//add(&v2, pos, left);
	add(&v2, v1,  direction);

	add(&v3, pos, direction);
	add(&v3, v3,  direction);

	add(&v4, pos, right);

	//add(&v5, pos, right);
	add(&v5, v4,  direction);

	vertices[offset+0]=float_to_vertex(c, up, make_float2(0,0));
	vertices[offset+1]=float_to_vertex(pos, up, make_float2(0,0));
	vertices[offset+2]=float_to_vertex(v1, up, make_float2(0,0));
	vertices[offset+3]=float_to_vertex(v2, up, make_float2(0,0));
	vertices[offset+4]=float_to_vertex(v3, up, make_float2(0,0));
	vertices[offset+6]=float_to_vertex(v4, up, make_float2(0,0));
	vertices[offset+5]=float_to_vertex(v5, up, make_float2(0,0));
	
}

///Generate leafs for the branches located in b array
/// for each branch n leafs is created
__global__ void generateL(SimpleVertex *vertices, DWORD* indices, branch * b, int n, float size)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
  
	
	
	branch current_b=b[idx];
	

	int offset_v=idx * n * VPL; 
	int offset_i=idx * n * IPL;


	
	float3 direction= make_float3(current_b.direction[0], current_b.direction[1], current_b.direction[2]);
	float3 position= make_float3(current_b.position[0], current_b.position[1], current_b.position[2]);
	normalize(&direction);


	float3 side;
	if(acosf(dot(direction, make_float3(0,1,0)))>0.1f)
		side=cross(direction, make_float3(0,1,0));
	else
		side=cross(direction, make_float3(1,0,0));
	normalize(&side);
	float3 up=cross(direction, side);
	normalize(&up);

	float step=current_b.length / n; //distance between leafs on the branch
	
	for(int i=0; i<n; i++)
	{
	 
		float3 d;
		float3 p;
		mul(&d, direction, step*i);
		add(&p, position, d);
		rotate(&up, up, direction, 2); /// 2 = random rotation

		vertices[offset_v+i*VPL]=float_to_vertex(p, up, make_float2(0,0));

		

		indices[offset_i+i*IPL]=offset_v+ i*VPL; 


	
		
		
		
	}


	
}





/************************************************************************/
/* Tree generation   */
//b = branches data
//s= numebr of branches
//w= number of m_walls per branch
//r = branch radius
/************************************************************************/
int generateTree(branch * b, int s, int w, float r,  SimpleVertex* vertices, DWORD * indices , ofstream &out)
{

	if(!InitCUDA()) {
		return 1;
	}


	for(int i=0; i<s; i++)
		out<<b[i].position[0]<<" "<<b[i].position[1]<<endl;


	int N = s * w * VPW;   //number of vertices
	int I = s * w * IPW;	 //number of indices
	size_t size = N * sizeof(SimpleVertex);
	size_t size_i = I * sizeof(DWORD);
	size_t size_b = s*sizeof(branch);
	
	SimpleVertex * d_v;  //device array with vertex data
	DWORD		 * d_i;  //device array with indices data
	branch	     * d_b;
	
	
	CUDA_SAFE_CALL(cudaMalloc((void **) &d_v, size));   // Allocate array on device	
	CUDA_SAFE_CALL(cudaMalloc((void **) &d_i, size_i));   // Allocate array on device
	CUDA_SAFE_CALL(cudaMalloc((void **) &d_b, size_b));   // Allocate array on device

	out<<(cudaGetErrorString(cudaGetLastError()));
	
	// Initialize host array and copy it to CUDA device
	
	CUDA_SAFE_CALL(cudaMemcpy(d_v, vertices, size, cudaMemcpyHostToDevice));
	CUDA_SAFE_CALL(cudaMemcpy( d_i, indices, size_i, cudaMemcpyHostToDevice));
	CUDA_SAFE_CALL(cudaMemcpy( d_b, b, size_b, cudaMemcpyHostToDevice));
	
	
	
	
	// Do calculation on device:
	int block_size = 2;
	int n_blocks = s/block_size + (s%block_size == 0 ? 0:1);
	// Print results
	out<<endl<<"before"<< N<<" "<<I<<" "<<w<<" "<<endl;
	generate <<< n_blocks, block_size >>> (d_v, d_i, d_b, w, r);
	CUDA_SAFE_CALL(cudaThreadSynchronize());
	

	
	// Retrieve result from device and store it in host array
	CUDA_SAFE_CALL(cudaMemcpy(vertices, d_v,  sizeof(SimpleVertex)*N, cudaMemcpyDeviceToHost));
	CUDA_SAFE_CALL(cudaMemcpy(indices, d_i,  sizeof(DWORD)*I, cudaMemcpyDeviceToHost));
	
	// Print results
	for (int i=0; i<N; i++) out<< i<<" "<<vertices[i].Pos.x<<" "<<vertices[i].Pos.y<<endl;;
	out<<(cudaGetErrorString(cudaGetLastError()));
	
	// Cleanup
	cudaFree(d_v);
	cudaFree(d_b);
	cudaFree(d_i);
	
	return 0;
}


int generateLeafs(branch * b, int s, int n, float leaf_size,  SimpleVertex* vertices, DWORD * indices , ofstream &out)
{

	if(!InitCUDA()) {
		return 1;
	}


	int N = s * n * VPL; 
	int I = s * n * IPL;
	size_t size = N * sizeof(SimpleVertex);
	size_t size_i = I * sizeof(DWORD);
	size_t size_b = s*sizeof(branch);
	
	SimpleVertex * d_v;  //device array with vertex data
	DWORD		 * d_i;  //device array with indices data
	branch	     * d_b;
	
	
	cudaMalloc((void **) &d_v, size);   // Allocate array on device	
	cudaMalloc((void **) &d_i, size_i);   // Allocate array on device
	cudaMalloc((void **) &d_b, size_b);   // Allocate array on device

	out<<(cudaGetErrorString(cudaGetLastError()));
	out<<endl<<"LISCIE"<<endl;
	
	// Initialize host array and copy it to CUDA device
	
	cudaMemcpy(d_v, vertices, size, cudaMemcpyHostToDevice);
	cudaMemcpy( d_i, indices, size_i, cudaMemcpyHostToDevice);
	cudaMemcpy( d_b, b, size_b, cudaMemcpyHostToDevice);
	
	
	
	
	// Do calculation on device:
	int block_size = 2;
	int n_blocks = s/block_size + (s%block_size == 0 ? 0:1);
	generateL <<< n_blocks, block_size >>> (d_v, d_i, d_b, n, leaf_size);
	cudaThreadSynchronize();
	

	
	// Retrieve result from device and store it in host array
	cudaMemcpy(vertices, d_v,  sizeof(SimpleVertex)*N, cudaMemcpyDeviceToHost);
	cudaMemcpy(indices, d_i,  sizeof(DWORD)*I, cudaMemcpyDeviceToHost);
	
	// Print results
	for (int i=0; i<N; i++) out<< i<<" "<<vertices[i].Pos.x<<" "<<vertices[i].Pos.y<<" "<<vertices[i].Pos.z<<endl;
	//out<<(cudaGetErrorString(cudaGetLastError()));
	
	// Cleanup
	cudaFree(d_v);
	cudaFree(d_b);
	cudaFree(d_i);
	
	return 0;
}



